//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry                         : env klmperry

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"

// base
cd  $perrydatas
use perry_base.dta, clear
gen     male  = 1 if C4 == 1
replace male  = 0 if C4 == 2 
gen     treat = 1 if C2 == 2
replace treat = 0 if C2 == 1
gen perry = 1
keep    id treat sb_iq_5 sb_iq_7 sb_iq_10 male perry
rename sb_iq_5  iq0 
rename sb_iq_7  iq2 
rename sb_iq_10 iq5
gen iq15 =0
tempfile perry
save "`perry'", replace

// abc
cd $abcjpeanalysis
use append-abccare_iv.dta, clear
keep if program == "abc"
egen iq0  = rowmean(sb4y sb5y)
egen iq5 = rowmean(sb5y iq12y iq8y iq15y)
egen iq2 = rowmean(iq6y6m iq7y iq6y)
rename iq21y iq15
keep id treat iq0 iq2 iq5 iq15 male
gen abc = 1
tempfile abc
save "`abc'", replace

// ihdp
cd $dataihdp

use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen     sample1 = e(sample)
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

global inputs_std   cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg sample $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs ppvtstd wasifsiq

// merge in ages 5 and 8
merge 1:1 ihdp using iqchildhoodihdp
keep  if _merge == 3
drop _merge

// impute iq
replace kidppvt5 = wippsif5 if kidppvt5 ==.
replace wippsif5 = kidppvt5 if wippsif5 ==.

replace ppvtstd8 = fsiq8    if ppvtstd8 ==.
replace fsiq8    = ppvtstd8 if fsiq8    ==.

replace ppvtstd  = wasifsiq if ppvtstd  ==.
replace wasifsiq = ppvtstd  if wasifsiq ==.

reg fsiq8 ppvtstd8 kidppvt5 wippsif5 wasifsiq ppvtstd
gen sample2 = e(sample) 

// outcomes
merge 1:1 ihdp using adultihdp
keep if _merge == 3
drop    _merge

// construct "positive" outcomes
replace idle18y         = 1 - idle18y
replace sch_eversped18y = 1 - sch_eversped18y
replace sch_math18y     = 1 - sch_math18y
replace sch_read18y     = 1 - sch_read18y
replace sch_thrp18y     = 1 - sch_thrp18y
replace teen18y         = 1 - teen18y

gen     smoke18 = .
replace smoke18 = 1 if cigs18y >  0  & cigs18y != .
replace smoke18 = 0 if cigs18y == 0 & cigs18y != .
replace smoke18 = 1 - smoke18

gen     absence18 = .
replace absence18 = 0 if sch_dabs18y < 5  & sch_dabs18y != .
replace absence18 = 1 if sch_dabs18y >= 5 & sch_dabs18y != .
replace absence18 = 1 - absence18 

// marking missing
global ed18  sch_eversped18y sch_math18y sch_read18y sch_test18y
global beh18 smoke18 idle18y sch_thrp18y teen18y
egen missing     = rowmiss($ed18 $beh18)

// average
egen avg_outcome_educ18 = rowmean($ed18)        if missing == 0
egen avg_outcome_beh18  = rowmean($beh18)       if missing == 0
egen avg_outcome18      = rowmean(avg_outcome_educ18 avg_outcome_beh18) if missing == 0

// define sample
reg     avg_outcome18 if sample2 == 1
replace sample2 = 0   if e(sample) == 0
replace sample2 = 0   if   sample1 == 0

egen iq15 = rowmean(ppvtstd wasifsiq)
egen iq2  = rowmean(kidppvt5 wippsif5)
egen iq5  = rowmean(fsiq8 ppvtstd8)

foreach var of varlist iq2 iq5 iq15 {
	replace `var' =. if sample2 ==0
}
rename ihdp id
rename tg treat
rename iqcage iq0
gen    ihdp  = 1  
gen    ihdpt = 1 if ihdp == 1 & twin == 1
gen    ihdps = 1 if ihdp == 1 & twin == 0
append using "`perry'"
append using "`abc'"
gen control = 1 - treat
replace male = 1 if sex == 1 & ihdp == 1
replace male = 0 if sex == 0 & ihdp == 1 

// strata
global perry_strata
global abc_strata
global ihdpt_strata strata(bwg site treat)
global ihdps_strata strata(bwg site treat)

foreach prog in abc perry ihdpt ihdps { 
	gen `prog'm = 1 if `prog' == 1 & male == 1
	gen `prog'f = 1 if `prog' == 1 & male == 0
}

// treatment effects
// programs
foreach prog in perry perrym perryf abc abcm abcf ihdps ihdpsm ihdpsf ihdpt ihdptm ihdptf {
	foreach var of varlist iq0 iq2 iq5 iq15 {
	bootstrap, ${`prog'_strata} reps($bootstraps): reg `var' treat if `prog' == 1
	matrix b = e(b)
	matrix V = e(V)
	matrix se`prog'_`var' = sqrt(V[1,1])
	matrix md`prog'_`var' = b[1,1]
	matrix  t`prog'_`var' = abs(md`prog'_`var'[1,1]/se`prog'_`var'[1,1])
	matrix df`prog'_`var' = e(N) - e(rank)
	matrix  p`prog'_`var' = 2*(1 - normal(t`prog'_`var'[1,1]))
	matrix  n`prog'_`var' = e(N)
	
	reg `var' control treat if `prog' == 1, nocons
	matrix b = e(b)
	matrix bc_`var' = b[1,1]
	matrix bt_`var' = b[1,2]
	
	matrix `prog'_`var' = [bc_`var',bt_`var',md`prog'_`var',p`prog'_`var',n`prog'_`var']
	}
	matrix `prog' = [`prog'_iq0 \ `prog'_iq2 \ `prog'_iq5 \ `prog'_iq15]
	matrix colnames `prog' = bc_`prog' bt_`prog' md_`prog' p_`prog' n_`prog' 
}

// graph
matrix all = [perry,perrym,perryf,abc,abcm,abcf,ihdps,ihdpsm,ihdpsf,ihdpt,ihdptm,ihdptf]
clear
svmat all, names(col)
foreach var of varlist *perry {
	replace `var' =. if _n == 4
}
// labels
foreach prog in perry perrym perryf abc abcm abcf ihdps ihdpsm ihdpsf ihdpt ihdptm ihdptf {
	tostring md_`prog', gen(md_`prog'str) force format(%15.1fc)
	replace  md_`prog'str = "{&Delta} = " + md_`prog'str
}

gen ppt = 105
gen nperry = _n*4 - 3
gen nabc   = nperry + 1
gen nihdp  = nabc   + 1
replace nabc  =  nabc - 1 if _n == 4
replace nihdp = nihdp - 1 if _n == 4

cd $output
#delimit
twoway     (bar     bc_perry  nperry if _n <= 3, color(gs10) barwidth(.8)              yline(105, lpattern(solid) lcolor(gs14)))
	   (bar     bt_perry  nperry if _n <= 3, fcolor(none) lcolor(gs0) barwidth(.8) yline(65, lpattern(solid) lcolor(gs14)))
	   (scatter ppt  nperry if _n <= 3, msymbol(none) mlabel(md_perrystr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_perry    nperry    if  p_perry < 0.05 & _n <= 3                     , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_perry    nperry    if  p_perry >= .05 & p_perry    <= 0.10 & _n <= 3, msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_abc  nabc, color(gs10) barwidth(.8)             )
	   (bar     bt_abc  nabc, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nabc, msymbol(none) mlabel(md_abcstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_abc    nabc    if  p_abc < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_abc    nabc    if  p_abc >= .05 & p_abc    <= 0.10 , msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_ihdps  nihdp, color(gs10) barwidth(.8)             )
	   (bar     bt_ihdps  nihdp, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nihdp, msymbol(none) mlabel(md_ihdpsstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_ihdps    nihdp    if  p_ihdps < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_ihdps    nihdp    if  p_ihdps >= .05 & p_ihdps    <= 0.10  , msize(small)  mcolor(gs0) msymbol(square))
	  
	   
	   ,
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(65[10]105, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Perry" " " " " "' 2  `" "ABC" " " "{bf:0}" "'  3  `" "IHDP" "(Singletons)" " " "'
			 5 `" "Perry" " " " " "' 6  `" "ABC" " " "{bf:2}" "'  7  `" "IHDP" "(Singletons)" " " "'
			 9 `" "Perry" " " " " "' 10 `" "ABC" " " "{bf:5}" "'  11 `" "IHDP" "(Singletons)" " " "'
			 13 `" "ABC" " " " " "' 13.5 `" " " " " "{bf:15}" "' 14  `" "IHDP" "(Singletons)" " " "'
		  
		  ,  labsize(small) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  xtitle("{bf: Years After Program Ends}")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export cskills_seriessingletons_all.eps, replace

#delimit
twoway     (bar     bc_perry  nperry if _n <= 3, color(gs10) barwidth(.8)              yline(105, lpattern(solid) lcolor(gs14)))
	   (bar     bt_perry  nperry if _n <= 3, fcolor(none) lcolor(gs0) barwidth(.8) yline(65, lpattern(solid) lcolor(gs14)))
	   (scatter ppt  nperry if _n <= 3, msymbol(none) mlabel(md_perrystr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_perry    nperry    if  p_perry < 0.05 & _n <= 3                     , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_perry    nperry    if  p_perry >= .05 & p_perry    <= 0.10 & _n <= 3, msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_abc  nabc, color(gs10) barwidth(.8)             )
	   (bar     bt_abc  nabc, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nabc, msymbol(none) mlabel(md_abcstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_abc    nabc    if  p_abc < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_abc    nabc    if  p_abc >= .05 & p_abc    <= 0.10 , msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_ihdpt  nihdp, color(gs10) barwidth(.8)             )
	   (bar     bt_ihdpt  nihdp, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nihdp, msymbol(none) mlabel(md_ihdptstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_ihdpt    nihdp    if  p_ihdpt < 0.05                       , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_ihdpt    nihdp    if  p_ihdpt >= .05 & p_ihdpt    <= 0.10  , msize(small)  mcolor(gs0) msymbol(square))
	  
	   
	   ,
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(65[10]105, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Perry" " " " " "' 2  `" "ABC" " " "{bf:0}" "'  3  `" "IHDP" "(Twins)" " " "'
			 5 `" "Perry" " " " " "' 6  `" "ABC" " " "{bf:2}" "'  7  `" "IHDP" "(Twins)" " " "'
			 9 `" "Perry" " " " " "' 10 `" "ABC" " " "{bf:5}" "'  11 `" "IHDP" "(Twins)" " " "'
			 13 `" "ABC" " " " " "' 13.5 `" " " " " "{bf:15}" "'  14 `" "IHDP" "(Twins)" " " "'
		  
		  ,  labsize(small) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  xtitle("{bf: Years After Program Ends}")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export cskills_seriestwins_all.eps, replace

#delimit
twoway     (bar     bc_perrym  nperry if _n <= 3, color(gs10) barwidth(.8)              yline(105, lpattern(solid) lcolor(gs14)))
	   (bar     bt_perrym  nperry if _n <= 3, fcolor(none) lcolor(gs0) barwidth(.8) yline(65, lpattern(solid) lcolor(gs14)))
	   (scatter ppt  nperry if _n <= 3, msymbol(none) mlabel(md_perrymstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_perrym    nperry    if  p_perrym < 0.05 & _n <= 3                     , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_perrym    nperry    if  p_perrym >= .05 & p_perrym    <= 0.10 & _n <= 3, msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_abcm  nabc, color(gs10) barwidth(.8)             )
	   (bar     bt_abcm  nabc, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nabc, msymbol(none) mlabel(md_abcmstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_abcm    nabc    if  p_abcm < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_abcm    nabc    if  p_abcm >= .05 & p_abcm    <= 0.10 , msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_ihdpsm  nihdp, color(gs10) barwidth(.8)             )
	   (bar     bt_ihdpsm  nihdp, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nihdp, msymbol(none) mlabel(md_ihdpsmstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_ihdpsm    nihdp    if  p_ihdpsm < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_ihdpsm    nihdp    if  p_ihdpsm >= .05 & p_ihdpsm    <= 0.10  , msize(small)  mcolor(gs0) msymbol(square))
	  
	   
	   ,
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(65[10]105, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Perry" " " " " "' 2  `" "ABC" " " "{bf:0}" "'  3  `" "IHDP" "(Singletons)" " " "'
			 5 `" "Perry" " " " " "' 6  `" "ABC" " " "{bf:2}" "'  7  `" "IHDP" "(Singletons)" " " "'
			 9 `" "Perry" " " " " "' 10 `" "ABC" " " "{bf:5}" "'  11 `" "IHDP" "(Singletons)" " " "'
			 13 `" "ABC" " " " " "' 13.5 `" " " " " "{bf:15}" "' 14 `" "IHDP" "(Singletons)" " " "'
		  
		  ,  labsize(small) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  xtitle("{bf: Years After Program Ends}")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export cskills_seriesmsingletons_all.eps, replace

#delimit
twoway     (bar     bc_perryf  nperry if _n <= 3, color(gs10) barwidth(.8)              yline(105, lpattern(solid) lcolor(gs14)))
	   (bar     bt_perryf  nperry if _n <= 3, fcolor(none) lcolor(gs0) barwidth(.8) yline(65, lpattern(solid) lcolor(gs14)))
	   (scatter ppt  nperry if _n <= 3, msymbol(none) mlabel(md_perryfstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_perryf    nperry    if  p_perryf < 0.05 & _n <= 3                     , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_perryf    nperry    if  p_perryf >= .05 & p_perryf    <= 0.10 & _n <= 3, msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_abcf  nabc, color(gs10) barwidth(.8)             )
	   (bar     bt_abcf  nabc, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nabc, msymbol(none) mlabel(md_abcfstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_abcf    nabc    if  p_abcf < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_abcf    nabc    if  p_abcf >= .05 & p_abcf    <= 0.10 , msize(small)  mcolor(gs0) msymbol(square))
	  
	   (bar     bc_ihdpsf  nihdp, color(gs10) barwidth(.8)             )
	   (bar     bt_ihdpsf  nihdp, fcolor(none) lcolor(gs0) barwidth(.8))
	   (scatter ppt  nihdp, msymbol(none) mlabel(md_ihdpsfstr)    mlabposition(12) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter bt_ihdpsf    nihdp    if  p_ihdpsf < 0.05                      , msize(small)  mcolor(gs0) msymbol(circle))
	   (scatter bt_ihdpsf    nihdp    if  p_ihdpsf >= .05 & p_ihdpsf    <= 0.10  , msize(small)  mcolor(gs0) msymbol(square))
	  
	   
	   ,
		  legend(label(1 "Control") label(2 "Treatment")
			 label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
			cols(2) rows(2) order(1 2 4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(65[10]105, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(1 `" "Perry" " " " " "' 2  `" "ABC" " " "{bf:0}" "'  3  `" "IHDP" "(Singletons)" " " "'
			 5 `" "Perry" " " " " "' 6  `" "ABC" " " "{bf:2}" "'  7  `" "IHDP" "(Singletons)" " " "'
			 9 `" "Perry" " " " " "' 10 `" "ABC" " " "{bf:5}" "'  11 `" "IHDP" "(Singletons)" " " "'
			 13 `" "ABC" " " " " "' 13.5 `" " " " " "{bf:15}" "' 14 `" "IHDP" "(Singletons)" " " "'
		  
		  ,  labsize(small) grid glcolor(gs14) glpattern(solid))
		  ytitle("Unadjusted Mean", size(medsmall))
		  xtitle("{bf: Years After Program Ends}")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export cskills_seriesfsingletons_all.eps, replace
